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Abstract. Wo propose a hierarchy of multi-level kinetic Monte Carlo methods for sampling high- 
dimensional, stochastic lattice particle dynamics with complex interactions. The method is based on 
the efficient coupling of different spatial resolution levels, taking advantage of the low sampling cost 
in a coarse space and by developing local reconstruction strategies from coarse-grained dynamics. 
Microscopic reconstruction corrects possibly significant errors introduced through coarse-graining, 
leading to the controlled-error approximation of the sampled stochastic process. In this manner, the 
proposed multi-level algorithm overcomes known shortcomings of coarse-graining of particle systems 
with complex interactions such as combined long and short-range particle interactions and/or com- 
plex lattice geometries. Specifically, we provide error analysis for the approximation of long-time 
stationary dynamics in terms of relative entropy and prove that information loss in the multi-level 
methods is growing linearly in time, which in turn implies that an appropriate observable in the 
stationary regime is the information loss of the path measures per unit time. We show that this 
observable can be either estimated a priori, or it can be tracked computationally a posteriori in the 
course of a simulation. The stationary regime is of critical importance to molecular simulations as 
it is relevant to long-time sampling, obtaining phase diagrams and in studying metastability prop- 
erties of high-dimensional complex systems. Finally, the multi-level nature of the method provides 
flexibility in combining rejection-free and null-event implementations, generating a hierarchy of algo- 
rithms with an adjustable number of rejections that includes well-known rejection-free and null-event 
algorithms. 
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1. Introduction. One of the widely used computational methods at atomistic 
scales simulating stochastic particle systems is the continuous time or kinetic Monte 
Carlo (kMC) method. The first implementations of rejection-free kinetic Monte Carlo 
methods in molecular simulations go back to the stochastic simulation algorithm 
(SSA) of Gillespie for well- mixed systems, [TT], and the n-fold method or the BKL 
method of Bortz, Kalos and Lebowitz, [3J, for spatially distributed Ising-type systems. 
Traditionally kMC algorithms are serial, explicit time-stepping methods, limiting their 
applicability, due to the high computational cost per event (time-step), to moderate 
size systems. The principal part of the cost consists of searching for an event and 
updating the reaction rates. In the last decade research works have been focusing on 
developing sophisticated search and update techniques to reduce this computational 
cost, [5l 1301 [33] . For example, search algorithms have been proposed using binary tree 
pointers [TO] [31] in order to obtain O (log N) complexity in the size N of the simu- 
lated system. The BKL method is designed to reduce the cost of the searching step by 
lumping the transition states into classes of the same probability. All these techniques 
even though they reduce the computational cost per event are highly demanding in 
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computer storage and implementation overhead. An alternative to the rejection-free 
methods is based on the uniformization of the simulated continuous Markov chain. 
The implementation leads to null-event algorithms, which at each step require calcu- 
lation of the rate for only one, arbitrarily chosen, event that is accepted or rejected 
with a probability provided by the uniformization of the process, [5]. Although null- 
event methods reduce significantly the computational cost per Monte Carlo step they 
can be highly inefficient when the acceptance probability of an event becomes small, 
resulting to exceedingly small time-steps. 

In this work we propose a partially rejection-free kinetic Monte Carlo method, 
the multilevel kinetic coarse-grained Monte Carlo (ML-KMC), for sampling high- 
dimensional lattice systems with complex interactions and/or lattice geometries. Our 
primary interest in simulating extended systems in which the system size N is large 
and the behavior is governed by the infinite volume or thermodynamic limit (i.e., 
regimes where N — > oo). We distinguish between long-range interaction potentials, 
whose range L ~ N 1 ^ is comparable to the system dimensions N 1 ^, where d is the di- 
mension of the lattice, and short-range potentials which have a fixed interaction range 
S independent of the system size N. Typically in the lattice-gas simulations the long- 
range potentials are Coulomb or Lenard- Jones and short-range ones are represented 
by nearest-neighbor or next nearest-neighbor potentials. Interactions that combine 
long and short-range potentials are difficult to be handled efficiently by existing BKL- 
type algorithms since the number of classes grow exponentially with the length of the 
interaction range and thus making an implementation intractable, [5|. A possible 
remedy is to simulate a system with compressed (coarse-grained) potentials that have 
a shorter interaction range. We have demonstrated that smooth long-range potentials 
can be coarse-grained with controlled errors leading to highly efficient coarse-grained 
Monte Carlo (CGMC) methods, [HI [HI M\ ED] . On the other hand a simple coarse- 
graining of singular, long-range as well as short-range potentials does not lead to an 
accurate approximation of the renormalized Hamiltonian and results to loss of im- 
portant microscopic information, due to the presence of strong short-ranged spatial 
correlations. As demonstrated in [21 [ [2] approximations of coarse-grained Hamilto- 
nians accounting for short-range or singular potentials necessarily involve multi-body 
interaction terms whose implementation can become quickly computationally expen- 
sive. The pitfalls of coarse-graining are well-known also from molecular simulations 
of polymeric systems where they produce spurious phase changes or incorrectly pre- 
dict existing phase transitions, [17] • Similarly, coarse-graining is expected to be 
challenging in systems with complex lattice geometries, e.g. [32], that can also induce 
complicated spatial correlations. 

In this paper we demonstrate the use of the proposed ML-KMC method in the 
efficient and accurate simulation of such complex systems where coarse-graining either 
fails or it is computationally very expensive. We focus on the important example of 
systems where competing short- and long- ranged interaction forces are present and 
lead to complicated phase diagrams and pattern formation in various physical and 
chemical systems, [2Tfl [5] \T7\. 

2. Overview of the proposed method. The key ingredient of the proposed 
method is the multi-level sampling of the evolution process based on the knowledge 
of an even less accurate coarse-grained, meso/macroscopic dynamics. The present 
study is an extension of j!2j from the equilibrium to dynamical sampling sharing the 
same principle of efficient coupling of different resolution levels. In [T2] we proposed a 
multi-level Coarse-Grainined Metropolis-Hastings algorithm appropriate for sampling 
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equilibrium properties of systems. Although the embedded Markov chain generated 
by the Metropolis algorithm converges to the correct equilibrium distribution it does 
not preserve physical dynamics, a fact that motivated the present work. We present 
the method in a general framework to demonstrate its applicability to on- and off- 
lattice systems, and present in detail the application to stochastic lattice systems with 
short- and long-range interactions. 

The dynamics of stochastic systems on a countable configuration space E are 
determined by a continuous time Markov process ({tr t } t >o, £), with the infinitesimal 
generator £ : L°°(E) ->■ L°°(E) defined by the rates c(cr, a'), a, a' e E: 

ww = c ^ ct ') - ^ ct )) > ( 2J ) 

for every observable defined as any cf> S L°°(E). More specifically, in kMC we typically 
compute expected values of such observables, that is quantities such as 

u((, t) := E< IfW)} = E fW)PW t; > (2-2) 

conditioned on the initial data <jq — (. On the other hand, the evolution of the 
entire system at any time t is described by the transition probabilities P{a,t;X) := 
P (at = a | Co — C) where C S E is any initial configuration. The transition probabili- 
ties satisfy the Forward Kolmogorov Equation (Master Equation), [9], 

d t PWt;Q= E cW,r)PW,t;0-cW<r') p W,t;0, (2-3) 

a' ,a'^a 

where P(<r, 0; C) = 5(a — Q and 5(a — () = 1 if a = C, and zero otherwise. By a 
straightforward calculation using (|2.3p we obtain that the observable (|2.2p satisfies 
the initial value problem 

d t u((,t) = £u((,t) , u(C,0) = /(C), (2.4) 

The numerical implementation of the evolution of the process is realized with the 
embedded Markov chain {X n }„>o, X n = a n s t with transition probabilities 

PWW) = C ^^, A(a)=5>(<^'), (2.5) 

where p(a, a') is the probability of a jump from the state a to a' . The residence times 
St a for which the system stays in the state a before a jump is distributed according 
to an exponential law with the parameter A(cr). 

Method. The proposed method generates an approximate process (Wt}t>o, C-) of 
the stochastic process ({<r t } t >o, £). It is based on projecting the microscopic space 
E into a coarse space E with less degrees of freedom and on the knowledge of a 
coarse rate function c(r],r}') which captures macroscopic information from c(cr, a'). 
We denote the coarse space variables r\ = Ta defined by a projection operator T : 
E — > E. For example, for stochastic lattice systems that we elaborate on in this work, 
approximate coarse rate functions are explicitly known from coarse graining (CG) 
techniques of |15[ll6j. In this work we analyze a two-level approach, i.e., coupling two 
configuration spaces E and E with different resolutions, while a multi-level extension 
can be considered analogously. The ML-KMC method consists of the following steps, 
a schematic description is demonstrated in Figure 
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Level 2: 



c(cr,(T / )~c(fT,a ,/ ) 



Ta=n 



Crf (<t'|t)',<t) 



Level 1: 



n 



c(n,n') 



Figure 2.1. Two-level decomposition, compressing with T and reconstructing with c r f, of the 
evolution process per event. 

i) Construct a computationally inexpensive approximating CG process on the coarse 
space £ described by a coarse generator £ with rates c(r], if). 

ii) Define the "reconstruction" rates c r f {<j'\rf , ct) constrained on the updated coarse 
state ?/ that are simple to simulate and such that 



The approximation and its error is quantified in Section [5] The overall procedure can 
be thought as a reconstruction in dynamics of stochastic processes from a CG process. 
Furthermore, this procedure generates stochastic processes that are controlled error 
approximations of the process ({ct}t>o, £), determined by the reconstruction rates 
c r f(a' 1 1]', a) and the level of coarsening. The function c r f (c'|t/, a) enriches the CG 
procedure by re-inserting details that were smoothed out by the coarsening proce- 
dure. Implementation. The multilevel nature of the method provides the flexibility of 

combining rejection- free and null-event implementation algorithms at each resolution 
and consists of: 

i) A rejection- free algorithm for sampling in the coarse space with 0(77,77') where 
the reduction of the computational cost compared to microscopic sampling is 
significant due to the compression of spatial scales and interactions range. The 
rejection-free algorithm selects the most probable coarse state 77' that the system 
will evolve. 

ii) A null-event algorithm for sampling at the fine space with c T {(a'\r]' , a), however 
with a low rejection rate due to the fact that 77' was chosen as the most probable 



This partially rejection- free implementation approach suggests a non-constant time 
step update in contrast to null-event methods where the time update is uniform for 
all system states. Non-constant time step updating algorithms have been proposed 
in [T] designing a class of kinetic Monte Carlo algorithms with an adaptive time step 
interpolating between BKL and null-event algorithm. The variable time step in these 
algorithms is based on the adaptively improved upper bounds of the exponential 
parameter controlling the time step distribution, while in our proposed method it is 
the result of the combination of BKL and rejection-free implementations on different 
resolution state spaces. However, one can further enhance the ML-KMC method with 
an adaptive time step kinetic Monte Carlo in the spirit of pQ, in view of the freedom 
on the choice of implementation techniques in each level. 



c(rj,r)')c r {(a'\r)' ,a) approximates c(a,a'). 



event. 



Analysis. We provide numerical analysis for (a) finite-time weak error estimates and 
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(b) long-time stationary dynamics for the proposed ML-KMC algorithm. The pri- 
mary challenge in the finite-time weak error estimates is to obtain bounds that are 
independent of the high-dimension of the interacting particle system, and this is ac- 
complished by focusing on suitable macroscopic observables. However, this technique 
leads to estimates that involve constants that grow exponentially in time, as is the 
case also in many of the classical numerical analysis estimates for stochastic differen- 
tial equations or partial differential equations that rely on Gronwall inequality-type 
arguments. In order to overcome this difficulty we develop a different approach that 
estimates the error in long-time behavior, i.e. in stationary regimes for reversible 
or irreversible processes, using the relative entropy of the path measure. In kMC 
as well as in molecular simulations in general, we are often primarily interested in 
long-time, stationary regimes, including (i) the behavior of the stationary measure 
and its sampling, as well as (ii) stationary dynamics, i.e. dynamics where the initial 
measure is the stationary distribution reached after long-time integration. The latter 
is an especially important regime describing dynamic transitions between metastable 
states in complex energy landscapes, while at this regime we construct the system's 
phase diagrams. 

The error analysis study we present reveals that the relevant quantity for assessing 
long-time simulations in the stationary regime is the information loss per unit time in 
the path space measure. In fact, the issues related to the error analysis in the station- 
ary dynamics regime is the primary novelty in our paper from a numerical analysis 
perspective and they are not restricted only to the multilevel kMC algorithms; they 
are ubiquitous for numerical approximations of reversible and irreversible stochastic 
dynamics such as Langevin stochastic differential equations (SDEs), thermostatted 
dynamics, dissipative particle dynamics (DPD) methods, etc. For such systems there 
is a wealth of approximating schemes, the simplest being in the time stepping, e.g., 
explicit, implicit, predictor-corrector, operator splitting, etc. We expect that our pro- 
posed entropy-based perspective could be used to assess such numerical schemes at 
the stationary dynamics regime in a quantifiable manner, in a variety of stochastic, 
extended as well as finite-dimensional systems. 

We begin this work in Section |3] presenting the method and continue with Sec- 
tion proposing an efficient implementation strategy leading to a partially rejection 
free method. Application of the method in stochastic lattice systems is presented in 
Section [4] for adsorption-desorption Arrhenius dynamics. In Section [5] we prove error 
estimates that provide a quantitative control of the approximating process. Section [7] 
establishes computational efficiency of the method over conventional sampling tech- 
niques. The performance of the method is tested in Section [5] using a benchmark 
model of Arrhenius dynamics with a competing short- and long-range interaction 
potential. 

3. Multilevel kinetic Monte Carlo. The ML-KMC method is a kinetic Monte 
Carlo method generating controlled- error approximate dynamics of the stochastic pro- 
cess ({at}, C)t>o, based on the decomposition of the rate function c(cr, a') into coarse 
and corresponding reconstructing terms, such that 

c(a, a') « c(a, a') = [] c«(%, r,'M?(4-iM, Vi-i) , 

i=0 

where r/o = a and rji = Ti_i7/i_i,i = 1 . . . , / are variables in a hierarchy of coarse 
spaces with the decreasing numbers of degrees of freedom. Coarsening consists of 
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projecting the microscopic space into a coarse space E with less degrees of freedom 
for which a coarse rate function c(t7, 77') is appropriately defined, where r\ G £ denotes 
the coarse space variables defined by a projection operator T : E —> E, Ter = 77. 

For the sake of simplicity we present the two-level method (ML-KMC) while 
every step in the study that follows can be easily adopted to the multilevel case. 
The construction of the ML-KMC, sketched in Figure [2j consists of two steps. In 
the first step we construct an approximating process on the coarse space E described 
by a generator C with rates 0(77,77'), extracting macroscopic information from the 
rates c(cr, a'). In the second step we construct rates c T i(a'\rj' ,a), simple to simulate, 
and such that 6(77, 7/)c r f {o~'\f]'i o~) approximates c(a,a') with an error quantified in 
Section [5j Following the above description we define a stochastic process with rates 

c(a, a') = 5(77, i]')cr{(o-'\i]' , a) , (3.1) 

that generates a corresponding continuous time Markov chain with transition proba- 
bilities 

p(a, a') = , A>) = £ c(<7, </) . (3.2) 

A ( Cr ) cr'GS 

The corresponding continuous time process ({dt}t>o, £) is defined by the generator 
C on E, 

Ctf>(a) = J] c(t, <0 (0(<O - 0(a)) , (3.3) 

cr'SE 

for every <f> G L°°(E). In the application we study in Section such coarse and recon- 
structing rates are explicitly defined, both for exact and for controlled error sampling 
approaches. In Section[S]we provide estimates that quantify the approximating errors, 
for finite and long-time regimes, with respect to the level of resolution and the inter- 
pretation of the rate function decomposition. Furthermore, controlled approximations 
such as (|3.1j) can have significant computational advantages, see Section [7] 

3.1. Partially rejection-free implementation of the ML-KMC method. 

Rejection- free methods are based on calculating (updating) rates c(a, a') for all a' G E 
at each Monte Carlo step and choosing (searching) an event that evolves the system's 
state based on the probabilities p(a,a'), see (|2.5p . In these methods every step pro- 
vides a successful event but the cost of implementation becomes formidable for large 
complex systems where the cost increases with the system size | E | and with the com- 
plexity of the transition rates. Uniformization provides a solution to this problem, 
suggesting a method that needs the calculation only of one transition probability 
at each step with the disadvantage of introducing rejections; a number of proposed 
events will not happen. Such methods are known as null-event methods, that generate 
a Markov jump process with the embedded Markov chain described by 

LiT'm 7>z (3 - 4) 

-j^-p(a,a'), if & 

where A* is a uniform upper bound of A(cr) in (|2.5[) . The time that the system stays at 
state er, r nu11 , is governed by an exponential distribution with the parameter A* > A(<r) 
for all a G E. As a result since 1/A* = E[t i1u11 ] < E[r CT ] = 1/A(«r) for all a G E, the null 
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event algorithm evolves the system with a smaller time step and needs more MC steps 
than a rejection free method. Despite this inefficiency, the significant reduction of the 
computational cost per MC step can be advantageous for high dimensional complex 
systems when compared to rejection free methods. The probability of rejecting a 
proposed state when the system is in the state a is 

<V) = 1 - ^ , (3-5) 

and it is controlled by A*, indicating that the tighter the upper bound is the less 
rejections are introduced. The combination of rejection- free and null-event techniques 
at the two levels of the ML-KMC method introduces a variety of sampling techniques. 
Here we propose, as the more efficient approach, the rejection-free method at the 
coarse level and the null-event algorithm at the microscopic level, balancing between 
an improved rejection rate of the null-event and the computational complexity of the 
rejection-free method. Given a £ S, ?y = Tcr the evolution of the system to a state 
a' £ E is achieved by the following: 

Method 1. ML-KMC 
Coarse level. Evolve to a new coarse state rf £ E with the probability 

P(v, v') = ~T7y~ > Hv) = ^2 cCn,v') ■ 

Microscopic level. Select randomly (uniformly) a' £ E under the constraint Tcr' = 
rj' , and accept it with the probability 

p T{ (a \rj,a) = — — — — , A rf (a) = max > c rf (a rj , a) , 

{(t':T.(t'— r)' ) 

or reject it with the probability 

{<t':T<t'=?/} 

For implementation specifics of this step we refer to Section [7] 
Time update. Update time by a random time step with an exponential law with 
parameter 

A*(cr)=A(r ? )A rf (cr). (3.6) 

With the following lemma we prove that the ML-KMC method provides correctly 
a partial uniformization of the rejection-free method. 
Lemma 3.1. For any a £ E we have 

X(a) < A*(ct). 



Proof. The statement follows directly from the simple calculation 

X(a) = £ c(a,a') = £ cfarf) ]T cy {a'\r,\ a) < A(^)A rf (cr) . 

o-'GS n'eT, {a>:Ta>= n >} 
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□ 

The real time updates, controlled by A*(cr), depend on the state of the system, 
while in Q3.4p the time step is uniform for all states, controlled by A*. The multilevel 
method has the rejection probability at the state <r, 

niulti/ n A(g) A(g) 

A(r?)A rf (a-) A*(cr) 

since 

cr'GS {cr':Tcr' =?) '} ' U 7 



CT , eE A(7?)A rf (<7) A(r?)A rf ((j) 



Next, we show that the rejection rate of the ML- KMC method can be controlled 
and depends on the approximation of the coarse-grained rates. Before stating the 
proposition we note that a process {at}t>o is defined as lumpable, ([!]), with respect 
to the coarsening procedure 77 = Tct when its rates satisfy 

£ c(<7, a') = 5(77, rf) . 

{ct':TV'=?7'} 

This relation implies uniform reconstruction rates c T {(a'\n' , a) = l/|{c' : Ter' = 77'} 
for all a 1 £ {a' : Tct' = 77'}, hence Y.{a':Ta'= v '} c rf ( a 'W> °) = 1> Kf( a ) = 1> 

A>) = J] c(a, a') = ]T J] ^ = E V) = X(V) , 

o-'ES r;'GS {<x':T<x'=r)'} jj'eE 

and 

multi/ \ -1 A(cr) 

A(j7)A rf (cr) 

On the other hand, the approximating process {crt}t>o hi ML-KMC is not necessarily 
lumpable, nevertheless its approximation error determines the rejection probability: 
Proposition 3.2. Let the coarse rates define an approximately lumpable process, 
that is 

Z(a,a')=c(r,,r ] , ) + 0(e), (3.8) 

{V:TCT'=r/'} 

uniformly in a,rj — Tcr, 77' for some e > 0. Then 

p™ m (a) = 0(e) . 



Proof. First, we denote by S I( = {77' : 0(77, 77') > 0}, that is all coarse configurations 
accessible in a single step from 77. Assumption (|3.8[) implies that 

5(77,7?') £ c T{ (a'\ V ',a) = c(r 1 ,r,') + 0(e), £ c r{ (a'\rj' ,a) ^ 1 + O(e) , 

{o-':Ter'=r/'} {a':Ta'=rf} 

(3.9) 
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hence A rf (cr) = max,/ J2{a':Ta'= v '} c rf WW, a ) = 1 + ®{ e )- Tnen 

A>) = J2 cfoO= E 5 M') E c rf (a'|r/,a) = ^(c(r7,r/) + 0( e )) 
= A(tj) + |E„|0(e) , 



and 

A>) _ A(r?) + |E„|0(e) = . ( 



A(77)A rf (a) A(»j)A rf (a) 

where in the last equality we have used that A r f(c) = 1 + 0(e) and A (77) ~ |E^| in 
view of the definition X(r]) = X^'es ^(77,77'). Therefore the rejection probability (I3.7P 
satisfies 

-multi 1 M 17 ) /nf,^ 



Remark 3.1. Other implementation strategies can be designed, for instance 
employing a rejection- free method at the second level and/or a null-event method 
at the first level. However, when A r f(er, 77') := X){er'-Ter'=jj'} Crf ( <7 'l 7 ''' <J ) ^ ne me thod 
can be implemented by a rejection-free algorithm at the microscopic level but the 
process generated will violate the Markovian property unless A r f(<r) = A r f(er, 77') for 
all 77' £ E. Another possible modification in the algorithm is performing a null-event at 
the coarse level combined with a null-event algorithm at the microscopic level. In this 
approach the need of a uniform normalization constant due to the null event nature in 
both steps is a disadvantage since then the time update is small. In fact the average 
time update will be inversely proportional to 

A* = N max{c(?7, 77')} max{c r f (cr \rj , a)} < A* , 

77,77' a' .a 

which indicates that the method will have higher rejection rate than a conventional 
null event algorithm. Even though the last two implementation techniques seem to 
have disadvantages, we expect that application in case studies could be effective, for 
example when A* ~ A* and/or the computational acceleration due to coarsening is 
significant. 

4. Applications to complex interacting particle systems. We demonstrate 
the use of the ML-KMC method in the efficient and accurate simulation of complex 
systems where coarse-graining either fails or it is computationally very expensive. In 
this paper we focus on the important example of systems with competing short- and 
long- ranged interaction; such systems typically have complicated phase diagrams, 
may exhibit pattern formation and arise in numerous physical and chemical systems, 
[29l |6j [27] . We discuss such a system next and use it as a demonstration example for 
the ML-KMC method in the following Sections. 

4.1. Microscopic dynamics. We consider an Ising-type system on a periodic 
(/-dimensional lattice A with N = n d lattice points. At each x G A we can define 
an order parameter a(x). For example, when taking values and 1, it can describe 
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vacant and occupied sites. The microscopic dynamics are described by a continuous 
time Markov chain with state space Sat = {0, 1} A . For a configuration a we denote by 
a x the configuration which differs from a by an order parameter flip at the site x. The 
configuration update a — ¥ a x occurs with a rate c(x, a), i.e., the order parameter at x 
changes over the time interval [t, t + At] with the probability c(x, a) At + o(At). The 
resulting stochastic process ({<7t}t>o, £■) is a continuous time Markov jump process 
with the generator defined in terms of the rate c(x, a) by (|2.1j) . In this work we present 
as an example the dynamics with the Arrhenius rate for the spin-flip (adsorption- 
dcsorption) mechanism 

c(x, a) = do (l — + doa(x) exp [ — f3U(x, a)] , (4.1) 

where U(x,a) = Hm{<j) — Hn(ct x ) is the interaction energy of a particle located at 
the lattice site x G A. Arrhenius laws are typically used in micro-kinetic modeling 
of physiochemical applications, see for instance [32] . The Hamiltonian function Hjy ■ 
Sat —¥ R defines the total energy of the configuration a G Ejy an d we consider pair- 
wise long- and short-range interactions 

H N {a) = ffW(a) + Jf«(tr) + £ h(x)a(x) , (4.2) 

where h = h{x) is an external field and 

= ^ J(x - y)a(x)a(y) , H « = ^ K(x - y)a(x)a(y) . 

We consider interaction potentials J, K : W l — > R depending on the distance of a; 
and y, where the distance \x — y\ does not have to be necessarily measured in the 
Euclidean norm. Furthermore, we assume: 

(Al) J{ x -y) = l i V^(^\x-y\) , x, 2 /GAandL>0 (4.3) 

V {1) £ C 1 (R) , V [l) (-r) = V {1) (r) , V® (r) = , for r > 1 
(A2) if G L 1 1 oc (R) and if (a; - y) ^ for \x - y| < 5 (4.4) 
and S is independent of the lattice size N. 

The parameters L and 5 then define a range of interactions, i.e., the number of 
particles interacting with a given particle at the site x. The parameter L can be 
equal to n, i.e., interactions with all particles. The scaling in (Al) ensures that in the 
infinite volume limit TV — > oo, J G L X (R). Note that the regularity condition imposed 
on V"' rules out singular potentials such as Coulomb, Lenard- Jones etc. However, 
such cases can be treated in our analysis by splitting the potential into a smooth part 
(with L = n) and the short-range singular part. The range parameter S <C n is fixed 
and independent of n. With this form of the Hamiltonian the energy difference can 
be expressed as 

U(x,a)= K(x-y)a(y)+ £ J(x-y)a(y)-h(x) = U^(x,a) + U^(x,a) , 

where the size of the support |l( s )| = S d = N (1) and with L ~ n we have = 
L d = O n {N), 
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4.2. Coarse grained dynamics. Systems with smooth long-range interactions 
are well approximated by coarse-graining techniques, [THl EH [H], and CGMC are 
reliable and highly efficient simulation methods with controlled error approximations. 
Furthermore, models where only short-range interactions appear are inexpensive to 
simulate with conventional methods as there exist algorithms with the complexity 
Oat(I) per time step, [25 . However, when both short and long-range interactions 
are present, the conventional methods become prohibitively expensive, while coarse- 
graining methods are either not easily applicable or very expensive to implement due 
to the necessity to incorporate in them multi-body interactions, [2"TI |2"]. 

In the series of papers [TH [15j [19] the authors initiated the development of math- 
ematical strategies for the coarse- graining (CG) in stochastic lattice dynamics. One 
constructs the coarse lattice Am by dividing A in M coarse cells, each of which con- 
tains Q = q d (micro-)cells. Each coarse cell is denoted by k £ Am- A typical 
choice for the coarse variable in the context of Ising-type models is the block-spin over 
each coarse cell Ck, 

V ■= \ v(k) = J2 a ( x ) : k e Am f ' 

I xEC k ) 

defining the coarse graining map T : — > Y*m, Tct = rj and the coarse state 
space Yjm = {0, 1, . . . , Q} Am . The coarse grained approximating process {r?t}t>o, 
[TBI [TBI |2"3"] , is defined by adsorption and desorption rates of a single particle in the 
coarse cell Ck 

c a {k,Ti) = do(Q - f](k)) , c d (k,T]) = d T}(k)exp[- pU(k,T))] , (4.6) 
where the CG interaction potential is given by 

u{k,ri)= ]T j(MM0 + ^(M)(»/(fc)-i)-fr(*)> 

where for the coarse cells k, I £ Am we define 

W x<£C k yeCi ^ W > xeC k y€C k ,y=£x 

as the interaction potential on the coarse space. 

4.3. Long-time behavior and the stationary measure. In many applica- 
tions of interacting particle systems the long-time behavior of (ergodic) evolution is 
characterized by the stationary (equilibrium) measure. If the rates c(x, a) satisfy the 
reversibility condition with respect to the measure iijq.p = Z^ l e~ l3HN ^' J \ also known 
as detailed balance, 

c(x, a)e- mN ^ = c(x, a x )e~ 0HN ^ , (4.8) 
then the jump dynamics leaves the Gibbs measure 

UnA**) = P N {da) (4.9) 



invariant, as is the case for (|4.ip . The factor Zn is the normalizing constant (parti- 
tion function). Furthermore, the product Bernoulli distribution P^{d(r), is the prior 
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distribution on A representing distribution of states in a non-interacting system. The 
total energy H^(a) of the system, at the configuration a — {u(x) : x G A}, is given 
by the Hamiltonian Hn (|4.2p . 

The coarse-grained process (|4.6I) satisfies detailed balance ensuring that the pro- 
cess, at least for long-range potentials J, has as its invariant measure an approximation 
of the coarse-grained Gibbs measure, |18j . 

fivAdl) = 4^e-^ (0) ^P M (dv) , (4.10) 



7(0) 



where 



H (0) {v) = -\ E E J(k,l)v(k)v(l)-y (0,0) E 1(0 (»7(0 - !) 
zsAm fceA M (eA M 

fc#! (4.11) 

+ E M^Mfc)- 

fceAM 

Applying the same coarse-graining formula (|4.7p to the short range potential K will 
introduce errors that are not well-controlled as N — > 00, [21 J. However, the multi- 
level technique provides an approach for constructing approximations that do not 
require higher-order (cluster) expansions of the short-range potentials developed in 
[21] . For use of the multi- level approach in the context of equilibrium sampling and 
corresponding error analysis we refer the reader to |12[ 113) . We also revisit this error 
analysis in Section [5l noting that we do not require the reversibility condition (|4.8j) 
for the simulated process in our results. 

4.4. ML-KMC method for Arrhenius dynamics. As a specific example we 
explain the ML-KMC method for sampling the microscopic process {cr t } 4 >o gener- 
ated by the Arrhenius rate for the adsorption-desorption mechanism. For the model 
considered the coarse space rate functions are explicitly given by the coarse graining 
technique in Section 14.21 The reconstruction rates that we present rely on approx- 
imate sampling with potential splitting where the long and short-range interactions 
are split to the first and second level respectively and no corrections to coarse-grained 
rates are applied. 

Approximate dynamics with potential splitting. In this variant only the sampling 
corresponding to the costly long-range interactions is performed at the first (coarse) 
level and the compression of the short-range interaction is avoided. The sampling of 
the short-range contributions is performed at the second (fine) level. The rates on the 
coarse space Sm 

c a (k, V ) = d (Q-v(k)) , c d (k,n) = dov(k)e- l30W ^\ 
where U®(k,T}) = E/eA M -/(MMO + J{k,k)(r)(k) - 1) - \h{k). Since a' = a x is a 
spin flip updating, the reconstruction rates in (|3.ip are explicitly defined and denoted 

by 

c&(*M) = i^, 4(* M = ^ e -/^W) , (4.12) 

where 

U^(x,a)= E K(x-y)a(y)~h(x),U^(x,a)= E J{x-y)<j{y)-\h{x) . 
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Note that c° f (x\k, 77) (and cj? f (x\k, 77)) are well-defined since T)(k) ^ Q (i](k) ^ 0) for all 
k E Am,! E Cfc when adsorption (desorption) process is selected in the cell k. With 
this choice of rates the ML-KMC method generates a Markov process ({<7 t } t >o,£) 
with the rate function defined by 

c(x,a) = c a {k,r\)c a Ti {x\k,r\) + c d (k, rj)c% (x\k, rj) 

= do(l - o-(a;)) + d cr(x)e-' , ^*> ,T 5 , (4.13) 

where we define 

U(x,a) = U {s) {x 1 a) + U {l) {k 1 T 1 ), xeC k ,r] = Ta. (4.14) 
In Appendix 1X1 we prove that c(x,a) satisfy the detailed balance condition with 

faA*r) = J-e-^ 5 "W P N (da) , (4.15) 
Zn 

and Zjv is the normalization constant corresponding to the Hamiltonian 

Hn{<t) = -|EE^- vM*Mv) - \ E E -W- «i/)M*Mi/) 

+ ^2h(x)a(x). (4.16) 

We define k(x) to be the coarse cell k £ Am such that a; E Cfe. 

5. Controlled-error approximations for complex systems. In this section 
we provide estimates that quantify the numerical error when approximating the con- 
tinuous time Markov process ({cr t } t >o, C) by {{a t }t>o, £) defined by (|4.13[) , with 
invariant stationary measures y,N,p{dcr), (|4.9p . and JiN,j3{dcr), (|4.15[) . respectively. We 
prove information loss estimates in long-time stationary regimes and weak error es- 
timates for suitably defined macroscopic observables in finite time intervals. Finally, 
as it is evident from the proofs, our results for stationary dynamics (processes) and 
weak error estimates, hold true for the ML-KMC approximation of general particle 
systems which are not necessarily reversible. 

5.1. Controlled approximations at long times. We analyze approximation 
properties in long-time, stationary regimes, including (a) the behavior of the station- 
ary measure and its sampling, as well as (b) the stationary process, i.e., dynamics 
where the initial measure is the stationary distribution reached after long-time inte- 
gration. The stationary dynamics present an especially important regime describing 
dynamic transitions between metastable states in complex energy landscapes, while at 
this regime we construct the system's phase diagrams, see the simulation of hysteresis 
in Figure [HH] 

The error analysis in the stationary regime is the primary novelty in our paper 
from a numerical analysis perspective and it is not restricted only to the ML-KMC 
algorithms; these questions are ubiquitous for numerical approximations of reversible 
and irreversible stochastic dynamics such as Langevin stochastic differential equa- 
tions, thermostated dynamics, dissipative particle dynamics methods, etc. We expect 
that our proposed entropy-based perspective could be used to assess such numerical 
schemes at the stationary dynamics regime in a quantifiable manner, in a variety of 
stochastic, extended as well as finite-dimensional systems. 
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(a) Estimates for the stationary measure. For reversible systems, the explicit knowl- 
edge of the invariant measures of the ML-KMC process allows us to compare them 
directly to the Gibbs states associated with reversible kMC. Error estimates are given 
in terms of the specific relative entropy 

n(fi\u) = N- 1 [ log {dti(a)/du(cr)}n(da) 

between the corresponding equilibrium Gibbs measures. The scaling factor iV -1 is 
related to the extensivity of the system, hence the proper error quantity that needs 
to be tracked is the loss of information per particle 72.(/zjv,^ |a*jv,^)- Relative entropy 
was used as measure of loss of information in coarse-graining in |T8l [20] , and as means 
for sensitivity analysis in the context of climate modeling problems, |28j . One of the 
results in [18] [22] concerns derivation of the loss of information per particle estimates 
on the coarse space for smooth long-range interactions J, (j4.3[) . In the following 
theorem we prove an analogous error estimate on the microscopic space taking into 
account both short and long-range interactions. 

Theorem 5.1. Let JiN,ff be the approximating measure of the microscopic equi- 
librium measure (J.n,/3 defined by l{4.15\ ) and |^.ff[ ), with Hamiltonian functions I{4-1()\ ) 
and h4-Zfy respectively, then the loss of information per particle is estimated by 



Before continuing with the proof of the theorem we state a necessary estimate in 
Lemma [531 that is proved in [18] . 

Lemma 5.2. Assume the interaction potential J satisfies (Al) in J^.ffp , then 
the coarse-grained interaction potential J, given by |^.7[ ) at the coarsening level q 
approximates for any x, y £ A and k, I 6 Am the potential J with the error 

\J{x-y)-J{k,l)\<2± sup ||W<V-J/')II <Cvjt, 

where the constant Cy is independent of q, L. 
Proof. [Theorem 15.1] 



i . z N 

— log ^— 

N & Z N 



{ d^ N ,p 
\dfJ,N,p 

4- — E~ 



/3(H N (a) - H N (a)) 



By the definition of the Hamiltonian and Lemma 15.21 we have the estimate 
±\H N (a) H N (a)\ = ±\H®(<r) - iT%)| < C^VV^U , 

where r\ = Ta and C is a positive constant independent of the system size N, the 
coarsening parameter q, and the short range potential K = K(x — y). Therefore 



1 1 ZN 1 1 TO 



and — 



exp{-/3(#jv(a) - H N (a))}] = O n (^|||VF«|| C 

O n (j3±\\VV®\\, 



/3(H N (o-) - H N (a)) 
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that concludes the proof. □ 

Theorem 15.11 proves that the potential splitting does not affect equilibrium prop- 
erties of the system, a fact that we indeed observe in the numerical experiment, see for 
example Figure 18.21 The error estimate is independent of the short-range interaction 
potential as was expected, since the approximation of the invariant measure results 
only from compressing the long-range interactions. 

(b) Approximating the stationary process. The analysis stems from attempting to 
understand the striking accuracy of ML-KMC in calculating phase diagrams, hystere- 
sis simulations, Figure I5TT1 and Figure 1531 as well as in dynamic transitions between 
metastable states, see Figure 18.41 and Figure 18.51 We assess the approximation of 
the kMC process {ct}te[o,T] by the ML-KMC {&t}te[o,T]i when the initial data are 
sampled from a stationary measure. We consider the relative entropy per particle 
formula in the time interval [0, T] 



K (V[ 0t T\\T>[o,T]) = J log 



dV 



[0,T] 



dV 



dV 



[0,T] 



[0,T] 



where 2?[ .t] (resp. T , [o,t\) i s the distribution of a process {o't}te[o,T\ (resp. {5t}tero,Tl) 
on the path space Q([0, T], Sjv), the space of right continuous with left limits Sat- 
valued functions defined on [0, T]. We note that the relative entropy measures the loss 
of information in the approximation of the kMC process {&t}te[o,T] by the ML-KMC 
{5t}te[u,T]- 

If the initial distribution of the process {ct}te[o,T] (resp. {5t}te[o,T]) is the station- 
ary measure p (resp. p), then the specific relative entropy simplifies to the following 
relation, [7J, 



ft(© [0 ,T]P[o,T]) =TH(V [0tT] \V [0 , T] )+K(p\p) , (5.1) 
where lZ(p\p) is the specific relative entropy between the stationary measures, and 



'K-(' I) [0,T]\'^ > [0,T}) 



— E u 

N 11 



A(a)-A>)-^A(<7)p(<x,<7')logC 



X(a)p(a, a') 
X{a)p(a, a') 



(5.2) 



is given in terms of the jump rates A, A and jump probabilities p, p of {o~t}te\o,T] 
and {&t}t£[o,T] respectively. Indeed, by Girsanov's formula, |24j . we obtain the corre- 
sponding Radon-Nikodym derivative 



dT> lo,T], , p(po) 

(P*)=^7-T ex P 



dV 



[0,T] 




Kp s )]ds-fYj>^p s )io g x Jp^^ 



X(a)p(a, p s ) 



dN s {p) 



(5.3) 

on any path {pt}te[o.T] in Q([0, T], Sat), where N s (p) is the number of jumps of the 
path p up to time s. Then for any continuous, bounded function (f> : S^r — > K 



cj>(p s ) dN s (p) 



= Ex> 


[L 







<j)(p s )X(p s ) ds 



TEMX], 



where Ex>[(j>(p s )} — J ( t>(Ps)dT>\ Q T ], the first equality results from the fact that N t — 
f„ X(p s ) ds is a (zero mean) martingale and the second equality follows because 
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{Pt}te[o.T] is a stationary process, i.e., Ex>[(f>(pt)} = E )l [(j)(p t )]. Thus we have 



K(©[o,r]P[o,T]) =N~ l E v 



log 



dV 



[0,T] 



[0,T] 



lo. 



+ N~ 1 E V 



f [X(p s )-X(p s )]ds-f ^p(a,p s )log ^l P J> Ps) dN s (p) 
Jo Jo aeTi X(a)p(a,p s ) 



A(a)-A(<7)-^A(aMa,a')log T 



X(a)p(cr, a') 



■k(m\P) 



X{a)p{a, a') 
= TH{V [0}T] \V [0iT] )+K(ti\fi) , 

which is the formula (|5.ip . 

Remark 5.1. Formula (|5 . 1 [) shows that in the stationary dynamics regime the 
information loss consists of two terms, one which scales as 0t(1) in T and is related to 
the stationary measures p and p and another one that captures the stationary dynam- 
ics and scales as 0t(T). Furthermore, we note that for the stationary process approx- 
imation the relevant quantity is the relative entropy per unit time H(X'[o ! t] \D[o,t]), 
(|5.2[) . On one hand, (15.11) implies that the loss of information increases linearly in 
time, while the stationary measure loss of information becomes irrelevant as T oo. 
The fact that as T grows in (|5.ip the term lZ(p\p) becomes unimportant is especially 
useful since p is typically not known explicitly (contrary to the case in (I4.15|l ). while in 
non- reversible systems (e.g., reaction-diffusion kMC) p is not known either, however, 
due to (|5.ip it is not necessary to calculate or estimate 1Z (p\p). 

We can use this information theory-based perspective to evaluate broad classes of 
numerical schemes for extended stochastic processes such as the ones arising in kMC, 
in the long-time, stationary process regime. However, here the following theorem 
provides the information loss estimates for approximating the microscopic process 
({°t}t£[o,T])£) with ({5f} te [ ,T]}£) defined by (|4.ip and f|4.13[) respectively. 

Theorem 5.3. [A priori estimates.] Let ({ci}te[o,T]> C) and ({at}te[o,T],£) be 
stationary processes with the initial distributions pn,/3 and pn,p respectively, then 
(a) For any fixed time T > 



K (£>[o,t]|£>[o,t] J = TH(T>[ 0tT ]\V[ 0iT ]) +TZ(p N ^\p N ^) 



(5.4) 



where 



n(V l0!T] \T) [0iT] ) = -E } 



A(a)-A(a)-^A(aMa,a')logC 



X(o-)p(a, a') 
X(a)p(a, a') 



(5-5) 

(b) For any N, coarsening parameter q < L and interaction potentials J(x — y), 
K(x — y) satisfying (Al), J^. 3fy, and (A2), J^.^[) ; respectively, 



H(p M \V [0 ,t] ) < PjC{K, V, 0) [[ W« ||i , 



(5.6) 



where || - [|i = \\-\\ L i is the L 1 norm on R andC{K,V,/3) = C^K^, \\V^ |U, 0) 
is a constant independent of N . 

Proof, (a) Relation (15.41) is a direct consequence of the earlier discussion. 
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(b) Recalling the definition of the Markov jump process for the microscopic and the ap- 
proximating process we have for all x G A, a G £jv that X(a)p(a, o~ x ) = c(x, cr), A(cr) = 
SaGA c ( x ' a ) an< ^ A(cr)p(2, a ) = ^(a;, cr'), A(cr) = X^eA^l 2 -' "^i with the rate functions 
c(x,cr) = d (l-cr(a;)) + docr(a;)e- /3C/ ( a: '' T ) andc(x,cr) = d (l -a(x)) + d (7(x)e- l3ij( - x '^ 
as defined in (|4.ip and (|4. 13|) respectively. Then according to formula (|5.2[) we have 



■H(D [0iT] \V [0t T}) = N- 



c(x, a) 



xeA 



c(x, a) 



We define A q ^ N (x, a) = U{x 1 a) — U(x, a). From the definition of U(x,a), (j4.14|) , 
it follows that A^foo-) = J7 (i) (^,o-) - U {l) (k, Ta), for all x G C fc . In view of this 
equality a straightforward application of Lemma 15721 states that there exists a constant 
c > such that the microscopic potential U(x,a) is approximated by U(x,o~) with 



Then 



\A q , N (x,a)\<cj-\\VV^\U for all cr G Sat, a; G A w 



(5.7) 



U(V [0tT] \V [0>T] ) = JV^E 



_ e -/9A giN (x,o-) 



^ d0Cr(s)e- /3C/(x ^ ) log ■ 



,(«)=! 



-0U(X,(T) 

-PU{x,tr) 



< Ar- 1 C(iT,y, y 8)E 



.a;G A 



AT 1 ! 



d a(x)e^ u ^A q ^(x,a) 



<2c^|c(ir,y, y 8)||vy«|| l! 

where C(K,V,P) = sup CT x exp{— /3U(x, a)}. 



Remark 5.2. /A posteriori error analysis] Reversing the roles of /j, and /2 in the 
formula (15.5[) we obtain an o posteriori calculation on the loss of information in (|5.4[) : 



n(f> [0 , T] \V [0iT] ) = ll A (A(ct) - A(cr)) - 53 A(a)p(cr, a') log 



AT 



A(cr)p(cr, cr') 



A(o>(cr, cr') 



(5.8) 



Indeed, viewing this function as an observable estimated on the approximating sta- 
tionary process, we note that it can be computed a posteriori in the course of an 
ML-KMG simulation by sampling from the stationary measure jl. We note that in 
[1"8] we derived and tested computationally a posteriori estimates for adaptive coarse- 
graining of extended systems, based on a similar relative entropy approach for sam- 
pling the stationary distributions. The a posteriori representation in (|5.8|) is general 
and applies to both reversible and irreversible processes and does not require the a 
priori estimates in Theorem 15.31 (b) (c) . The complexity of numerical calculation of 
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depends on the complexity of the studied model. For example, in the lattice 
systems we study here, the loss of information per unit time is sampled in the course 
of a ML-KMC simulation for T > 1 as 



1 n(v {0tT] \V l0 ,T)) ™H(V [0<T] \V [0iT] ) = ^E A [(A(o-)-A(o-))- ^ c(x,o-) log 



c(x, ct) 



c(x, ct) 

From the practical point of view we have to deal with the somewhat computationally 
costly summation over the lattice which, in principle, has the complexity of O(N). 

5.2. Weak error estimates in finite time. In this section we prove weak error 
estimates of the approximation for a class of macroscopic observable quantities defined 
below. The weak error is defined by e w = |E CTo [</>(<7t)] — E CTo [0(cr t )]| for an observable 
<fi on the microscopic configuration space Eat, where the expectation is defined for the 
path conditioned on the initial configuration cto. We provide a quantitative measure 
of the controllable approximation that depends on two features: (a) the coarsening 
level q and (b) the potential splitting. The explicit dependence on the strength of the 
short range interactions provides us a measure to control splitting of the interactions 
into short and long-range parts. 

Theorem 5.4. Let ({CT t } t > , C) be the Markov process generated by the con- 
ventional kinetic Monte Carlo method, and ({5t}oo, £>) the process generated by the 
ML-KMC method, both with initial condition cto- For any macroscopic observable, i.e. 
a function cf> £ L°°(T,]sf), such that when d x <j){a) :— 4>{<J X ) — </>(c) , 

||<9a;</>||oo < C < oo , where C is independent of N , (5-9) 

X 

the weak error satisfies, for < T < oo, 

|E CT0 [^(<7 T )] -E ffo [0(3F T )]| < C(K,V,P)C T ± (5.10) 

where Ct is a constant independent of the system size and C(K, V, /3) — KsJl, 
K s = \suv xa e-P u{s){x ^\ and J L = \ sup^ a e -^Hk( x ),Ta)y 

Some typical macroscopic observables satisfying (15.91) are the coverage, c(at) = 
7f 12x£A a t( x )i the spatial correlations /(ct; fc) = X^eA a ( x ) <J { x + h), and the 
Hamiltonian defined in (|4.2p . 

For the proof of the theorem we will need the following Lemma [5?5] that we prove 
in Appendix iBl see also [19]. We define u(t,ao) = E[</>(ctt)|ct4 = cto] and the function 
u(t, cto) solves the backward Kolmogorov equation, i.e., the final value problem 

d t u(t, a) + Cu(t, a) = , u(T,-) = <p, t<T (5.11) 

For all observables <f> satisfying (|5.9|) we can estimate d x u(t,a) — u(t,<r x ) — u(t,a) 
independently of N since we have the following estimate. 

Lemma 5.5. Let u(t, ct) a solution of i5.11]) where C is the infinitesimal generator 
£/(ct) = J2 X c ( x > a ){f{ aX ) ~ f( a )) defined by the rate function c(x, ct) given in j4-l\ ). 
Then for any t < T 

]rR U (i,-)iu <cv^R0iu (5.i2) 



We continue with the proof of Theorem 15.41 
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Proof. Using the martingale property we have for any smooth function v(t,ao) 
and the process {<J t }t>o with the generator C 

E CT >(T,a T )] = E ct >(0,ct )] + / E ao [{d s + C)v(s,a s ))}ds. 

Jo 

Therefore 

E^o [<t>(<rr)] - E* W)] = E* a NO, tr )] - E CTo [u(T, a T )] 

E CT0 [(3 S + £)u(s,a s )} ds 



, [Cu(s, a s ) - Cu(s, a s )] ds 



However, we can bound the last term by 



^ (c(x, tr s ) - c(x, <r a )) 9 K u(s, cr a ) 



,l£A 



ds . 



^ |9 x u(s,cr s )|<is 



We conclude by using Lemma 15.51 and noting that 

|c(x,? s )-c(x,a 5 )|<d oe -^ (!)(fc(:c) ' T5f)) 
where using (|5.7|) . 



e -/3a«(x,Ta„) _ e -^t/<"(x,5f s ) 



for some C > 0. □ 



< C e -^ (!) ( fe (*>' T? »|A ? ^(x,a)| < CJ L |-||W«|| 



6. Exact sampling of kMC dynamics. In addition to the approximate dy- 
namics discussed so far, the ML-KMC method can also generate the exact dynamics 
associated with the rates c(a, a') by appropriate choice of the coarse and reconstruc- 
tion rates, albeit at higher cost than the controlled-approximation dynamics c(cr, a'). 
More specifically, given c(a,a') and 0(77,77'), c r f (a' 1 77', a) can be selected such that 



c{n,r]')c r[ (a'\r]',a) = c(a,a') = c(a,a') 



(6.1) 



In this case the ML-KMC method generates exactly the same stochastic process with 
the direct kMC method achieving a perfect reconstruction. Relation (|6.1[) ensures 
that processes {5*}t>o and {<Tt}t>o have the same generator £ = £, which is suffi- 
cient to prove that the two processes arc identical, 26 . As a specific example we 
demonstrate exact sampling via the ML-KMC method for the microscopic process 
{o~t}t>o generated by the Arrhenius rate for the adsorption-desorption mechanism in 
Section HIT] For this model the coarse space rate functions are explicitly given by the 
coarse graining technique in Section 14.21 The reconstruction rates that we present 
rely on correcting the error introduced by coarsening at the null-event step. 



20 



E. Kalligiannaki, M. A. Katsoulakis and P. Plechac 



Indeed, the rates on the coarse space corresponding to the compressed inter- 
actions for U(x, a), (|4.5p . are given by 

c a {k, n) =d (Q- r)(k)) , c d (k, r)) = d o r,(k)e- 00 ^ , 

where 

u(k, v)= J2 + J( - k > 0M0 + fc ) + *)](»?(*) - !) - H k ) ■ 

The reconstruction rates are explicitly defined by 

&{x\k, „) = , rj) = , ( 6 .2) 

Q - r?(fc) r?(fc) 

where we can also compare them to the reconstruction of the approximate dynamics 
in (j4.12j) . This choice of rates ensures that the ML-KMC method generates the 
same process ({o~ t }t>o, £) since the two-level process has the rates c(x,o~) — c{x,a). 
Furthermore, the quantities U{x,o~) — U{k,rf) are better localized in the sense that 
they decay faster than U{x 1 cr), hence it is easier to compress through truncation, [2]. 
Finally, the detailed balance condition is satisfied with invariant measure fj,N,p{dcr)) 
for completeness we give the proof in Appendix \K\ 



Remark 6.1. In analogy to the path- wise decomposition (13.11) of the stochastic 
process, exact or controlled error equilibrium sampling has been achieved with mul- 
tilevel CGMC methods, [T3], based on an analogous decomposition of the sampling 
probability measure fj,(da), i.e., fi(da) = p,(dri)i'(dcr\rj), where v defines the recon- 
struction and ft is the measure on the coarse space. 

7. Acceleration and computational complexity. The purpose of this section 
is to compare the efficiency of a ML-KMC method with conventional methods. In the 
presence of long-range interactions sampling with a rejection-free algorithm is next 
to impossible due to the very high number of classes in BKL-type methods, see for 
instance Table 17.11 Hence we may inevitably be forced to use a highly inefficient 
null event algorithm. With the proposed ML-KMC approach, a rather crude CG 
of the long-range potential gives rise to much fewer classes, thus we can sample at 
the first level rejection-free using the BKL algorithm while the next level can be 
null-event. The ML-KMC algorithm is applied to the stochastic lattice model for 
Arrhenius dynamics in Section 01 where we consider the global search implementation 
as in the conventional stochastic simulation algorithm (SSA). 

Algorithm 1. Two-level ML-KMC 
Given a,r] = To- 
Coarse level. (Level I) 

Update, (a) Calculate transition rates c a {k, 77), Cd(k, rf), for all k € Am and 

Search. Obtain uniform random numbers U\,U2 G [0, 1). 

If \%[(r}) < Mi adsorb else desorb. Assume that adsorption is chosen, 
then find k 6 Am such that X^-iiv) < ^m( 7 ?) u 2 < Kiv)- 
Microscopic level. (Level II) 

Reconstruct. Pick uniformly a site x in the cell Cfc. 
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Accept/Reject. Select a uniform random number u G [0,1), and define 
c r f(a;|fc, 77) = Crf d (x\k, 77) according to the selected process in the coarse 
move, e.g. {Jljp or \6.2\) . If X T {(a)u < c r f (x\k, 77) accept and update 
the state at the site x. 
Time update. Update time from an exponential law with the parameter A*(er), with 

A*((7)=%)A rf (<7). 

Here A*(cr) = A(?y)A r f(o-) is defined according to the sampling strategy. Specifi- 
cally for the exact sampling of Section [SJ 

A(r?) = Y, Ml - V(k)) + d aV (k)e-^ k ^ , 

k 



A r f (c) = q max 




and for the approximating sampling of Section [4. 4[ 

Hv) = £ d a (q - V (k)) + dov(k)e^ mi)(k ^ , 

k 

Alf(a)=9maX {^W'^) e "^}' 

where U* = mm XtCr (U(x, a) — U{k,rf)) and II s * = mm X}Cr U^(x,a). 

The ML-KMC method provides an efficient balance between benefits and limita- 
tions of the conventional null-event and rejection-free methods, that we summarize 
in Table 17.11 This is achieved by (a) improving the computational cost of a conven- 
tional rejection-free method, see Table [821 and (b) increasing the successful events of 
a null-event method. An event is considered successful when it is accepted and the 
system evolves to a new state. The cost per event of a kMC algorithm can be divided 
in two categories. The search cost, the computational cost to choose an event, and 
the cost of updating the rates when an event is performed. In Table [7TT1 the updating 
cost is realized as the number of operations necessary to calculate energy differences 
appearing and the search cost as the length of the array from which the next event 
(site) is selected. We consider global search and update algorithms for the compar- 
ison here, but we mention that both the cost in the traditional rejection- free and 
ML-KMC algorithms can be improved with the use of a sophisticated search/update 
algorithm, for example, with binary tree methods, [5]. For the sake of comparison 
and completeness we describe the conventional sampling algorithms, SSA, BKL and 
null event, in Appendix [C] The last column of Table 17.11 reveals the acceleration of 
the method in generating successful events, for example, when the system is at the 
state a the rejection probability of a proposed event in the ML-KMC algorithm is 
given by (|3.7|) . 

Next we present another argument that reveals the fact that the proposed method 
improves the computational cost of kMC algorithms. The number of classes in a BKL 
algorithm is determined by the level sets of U(x,a), in fact we can write U{x, a) := 
U(k,Ta) + E(x,a). Thus the level sets of U(k,rj) are defined on a coarser lattice, 
hence U(x,a) has many more (q d more) level sets. For a potential decaying at the 
length L on the microscopic lattice, U(x,a) — ^2,J{x — y)o~(y) — 0(2 dL ) different 
values (and classes), while the coarse interaction potential decays at a distance L/q 
and U{k,Ta) = ^(MMO = 0{q dL / q ) = 0(2 dLlo s(«)/9) different values (and 
classes). Hence a BKL algorithm on the coarse space has, by the factor 2 dL ( 1 ~ i ° s ( q ^ q \ 
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less classes in the implementation of the BKL algorithm. Clearly when the range 
of interactions L is large, the number of classes grows exponentially with L and 
implementation of a microscopic BKL algorithm is not feasible. Therefore sampling 
with a null-event algorithm is unavoidable. 

Table 7.1 

Computational complexity and event rejection rate comparison for a single kMC step in one 
space dimension. 





Search 


Update 


Rejection rate 


Rejection free (SSA) 


0{N) 


0(L x L) 





Two-level ML-KMG (SSA) 


0{M) 


0(L/q X L/q) 


1-A(a)/A» 


Rejection free (BKL) 


0(2 L ) 


0(L x L) 





Two-level ML-KMC(BKL) 




0(L/q X L/q) 


1 -A(a)/A» 


Null - event 


0(1) 


O(L) 


1- A(<j)/A* 



The ML-KMC method achieves acceleration of the rejection free simulations up 
of order q when sampling for the same finite time interval. We also note that since 
a transition to a new event is based on a single spin-flip, the reconstruction step is 
performed locally, confined to a single coarse cell, a fact that improves further the 
computational cost of the method. The computational times (CPU) compared next 
are those needed for reaching the same real time T with the conventional SSA and ML- 
KMC method, where we consider that CPU time is proportional to the computational 
complexity of the algorithms given in Table 17.11 Let n be the number of MC steps 
necessary in a rejection- free method to reach real time T. The corresponding necessary 
MC steps in a ML-KMC method are m = n/E MJV „ [p^f(<r)], where p^(a) = 
1 — f>™j '(<t) is the acceptance probability of an event when the system is in the state 
a. For the search algorithm the cost ratio of the microscopic SSA and the ML-KMC 
method is 

_ CPU S!rc j-f r e C Nn _ r^muitii 

S ~ CPU s , multi Mm -^,,[Psuccj 

and for the update 

_ CPU M , re j-free L 2 U _ ~ [-multi! 

" ~ cw u , multi {L/ q y m ~ q ^ [Psucc 1 ■ 

8. Numerical experiments: an Ising-Curie- Weiss model. We consider a 
benchmark problem with competing short and long-range interactions that exhibits 
complex multi-phase behavior as captured in the phase diagrams in Figure 18.11 and 
Figure [8731 Exact solutions for the free energy in the thermodynamic limit, N — s> oo, 
are known for the one-dimensional and two dimensional models, [14] . The energy of 
the system at the configuration a = {a(x), x G Ajy} is defined by the Hamiltonian 

x \x—y\=l x y^x 

= H^ s \a)+H (l \a)+E(a). 

The interactions involved in H^ s '(a) are nearest-neighbor with the constant strength 
K, while H^'ia) represents the long-range interactions given by the potential J with 
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the range L — N and h is an external field. A closed form solution in the thermody- 
namic limit (N — > oo) for the total coverage cg(K, J, h) in the one-dimensional model 
was derived in [14] . 

c p (K, J, h) = \m p {\k, \j, \h-\j + (8-2) 

where Mp(K, J, h) is a solution (minimizcr) of the problem 

min [ ^m 2 - log(e K cosh(/i + JW) + (e 2K sinh 2 (/i + Jm)+r 2K ) 1/2 )l . 
m \2 J 

Depending on the system parameters (K, J, h) can be a multivalued function and 
phase transitions may occur. We are interested in sampling the dynamical behavior 
of the system in the bi-stable regimes as well as in constructing the phase diagram 
with respect to the external field h. 

The computational examples demonstrate both acceleration of simulations with 
ML-KMC and the improved accuracy of ML-KMC contrasted with the CGMC simu- 
lations. The reference solution is obtained by the fully resolved microscopic simulation 
performed by the traditional null-event kMC. In numerical implementations we tested 
the three methods discussed previously: 

(i) the direct null-event kMC, (Algorithm [4]) , 

(ii) the developed ML-KMC (Algorithm [1} sampling on the microscopic space, 

(iii) the null-event CGMC sampling on the coarse space only, i.e., with both short 

and long-range potentials coarse-grained and without corrections due to the 
reconstruction step in ML-KMC. 
The energy difference U(x,a) appearing in the transition rates (|4.1[) is given by 

N 

U(x,a)=K J2 v(y) + jY, a ^-' 1 - 

\m—y\=l y=l 

For the ML-KMC method we apply the potential splitting approach where the rates 
on the coarse space Em, at the first level of the method, are defined in (|4.6[) with the 
potential energy 

M , N , 

m\k,r,) = jY,ri(k)-~ = jY,<r(y)-~, 

fe=l ~ y=l 

and the reconstruction rates at the second level of the method are defined by (|4. 12[) 
with 

2-2/1 = 1 

To implement the null-event method we need a uniform upper bound of the rates 

Crf(a#,r?), mm, 

where Ta = rj and = \ min{0, K}\. Therefore the time step of the method, that 
is proportional to A*(<r) = X(rj) X r { (a) , clearly varies with the system state a since 
Mv) = j: k [do(q-v(k))+d r l (k)e-^ k ^]. 
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Note that in this example the coarse-grained Hamiltonian is exact, i.e., 

H®(rj) = if™((r), thus there is no approximation error due to coarse-graining the 
long-range potential. Therefore, while CGMC sampling is approximate, due to coarse- 
graining of both H^\rf) and H^[rj), the ML-KMC method samples the exact mi- 
croscopic process, i.e., c(x,a) = c{x,a) for all x G A, a G S^v. This allows us to 
quantify the effect of splitting the potential function into short and long-range parts. 
For example, Figure 18.71 shows that the potential splitting is not introducing errors, 
which verifies the theoretical estimate for the information loss of the equilibrium dis- 
tribution, Theorem 15.11 The effect of the splitting is apparent only in the average 
acceptance rate of the method where the strength of the short-range interactions K 
controls the rejection rate according to (|8.3[) . 

In order to test the effect of coarse-graining in the ML-KMC method we modify 
the long-range potential in the Hamiltonian (|8.1j) and consider finite-range interac- 
tions, with the range L < N , and with the long-range part of the Hamiltonian 

x \y—x\<L 

For this case the proposed ML-KMC method is approximate, however, it still reduces 
significantly the coarse-graining error of the direct CGMC sampling, since compressing 
of the short-range part is avoided, see Figure I57B1 and description below. 

In all simulations we consider the one-dimensional model with the coarsening 
parameter q = N in the CGMC and ML-KMC methods, that is the coarse space 
consists of one cell k = 1 and the coarse variable n is the total coverage 77 = To- = 

Stationary dynamics and equilibrium sampling. We demonstrate properties of the 
ML-KMC algorithm in the stationary regime by constructing (equilibrium) phase 
diagrams of the average coverage with respect to the external field h. We explore 
different regimes of the phase plane K-J . With the choice of the potential parameters 
K and J corresponding to bi-stable regimes we observe that the ML-KMC algorithm 
approximates properly the hysteresis behavior while the CGMC algorithm samples 
incorrect energy landscape and thus does not estimate the hysteresis behavior cor- 
rectly. The coarse-graining parameter is set to q = N both in ML-KMC and CGMC 
simulations. The ML-KMC method avoids compressing the short-range interactions 
that introduce large error in the CGMC simulations. Figure [57T1 depicts the hystere- 
sis behavior in the case when the long-range potential is of Curie- Weiss type, i.e., 
the interaction range is L = N, and hence it is coarse-grained exactly by block-spin 
coarse variables. However, coarse-graining the short-range, nearest-neighbor Ising, 
potential introduces an error which leads to a wrong prediction of hysteresis in the 
CGMC simulation. In Figure 15721 and 15751 we chose the interaction range L < N which 
also introduces coarse-graining error in the coarse-grained long-range potential. How- 
ever, the presented error analysis for the invariant measure suggests that this error is 
small and thus the ML-KMC sampling, unlike the CGMC simulations, are in a good 
agreement with estimates from the microscopic simulations. 

Transients and dynamical sampling. In the dynamical sampling we explore two quan- 
tities of interest: 

(a) The path-wise behavior of the coverage, defined by c(<7j) = jt YlxeA a t(%) for the 
direct microscopic sampling, c(at) — j? YlxeA &t{x) for the two-level sampling 
and c(rj t ) = j? J2keK M for the CGMC sampling. 
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Figure 8.1. Hysteresis simulation in a bi-stable regime. Potential parameters K = 3, J = 5, 
L = N , and the lattice-size N = 1024 and the coarsening parameter q = N . 

(b) The mean time to reach a transition from one equilibrium to another in the bi- 
stable regime, the exit time, t = E[T], T = inf{< > : c t > C}. The 
probability density functions (PDFs) p m , pa, p cg for the exit time estimated in 
the microscopic, the two-level ML-KMC and the CGMC methods respectively 
are monitored. Starting from an initial state with the coverage Co = in all 
methods we record the time r when the coverage exceeds the value C — 0.99. 
Estimating the observable r tests both a proper approximation of the energy land- 
scape as well as the correct time-scale in approximating dynamics. The simulation is 
set for the parameters K, J such that the system exhibits transition to an equilibrium 
which is, depending on the value of the external field h, stable or metastable (see Fig- 
ure 18. ip . We compare not only the expected (mean) values but also the probability 
density functions (PDFs) in order to demonstrate importance of error estimates in 
terms of the relative entropy. Probability density function was estimated from 10 4 
independent samples using the MATLAB estimator ksdensity with a normal kernel 
function. 

Figure [8T4l shows the comparison in the case of a single equilibrium state c pa 1.0 
(for the given value of h) and the long-range potential which is coarse-grained exactly, 
i.e., L = N. We observe a good agreement in all three methods, although the CGMC 
introduces a visible error due to coarse-graining of the short-range potential. 

An additional error is introduced by coarse-graining the long-range potential with 
L = 100 < N. Comparison of the exit time PDF in the case of a single equilibrium 
state c « 1.0 (for the given value of h) is depicted Figure 18.51 While ML-KMC 
simulations are in a good agreement with the microscopic simulation the CGMC 
algorithm introduces significant error for the estimated PDF. 

By adjusting the external field h the sampling is performed in the bi-stable regime 
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External field h 

Figure 8.2. Hysteresis simulation in a single phase regime. Potential parameters K = —5, 
J = 5, L = 20, the lattice size N = 256, and the coarsening parameter q = N . 

with two meta-stable equilibria. Coarse-graining both short and long-range potentials 
changes significantly the energy landscape and the CGMC algorithm cannot capture 
the transition within the simulation time-window. The exit-time probability distri- 
bution function is depicted in Figure 18.61 showing that the ML-KMC algorithm is 
capable of capturing the transition and approximate the exit-time PDF. The inset 
demonstrates that the CGMC simulation was unable to estimate the mean exit time 
as no transition occurred and the exit-time PDF is concentrated at the final time of 
the simulation window. This fact is further visualized in Figure I87H where evolution of 
the mean coverage is depicted. While the ML-KMC simulation results in a trajectory 
that approximates well the reference trajectory obtained from microscopic null-event 
kMC with a transition from c = to c w 1 equilibrium, the trajectory averaged in 
the CGMC simulation does not exhibit any transition in the simulation window. 

In Table [8TT1 we compare numerical results for the exit time and the correspond- 
ing computational times of the three algorithms for different values of the potential 
parameters. For the finite range L < N interactions, where coarse-graining error is 
present, we see that the ML-KMC method estimates are closer to the microscopic 
(conventional) method even when the CGMC method fails. Furthermore, we see 
a significant acceleration of the computational time both with the CGMC and the 
ML-KMC method. 

Appendix A. Detailed balance. 

A.l. Exact dynamics. The rate c(x,a), (|6.2p . satisfies the detailed balance 
condition with HN,p(d<j), (|4.9|) . i.e., 
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Figure 8.3. Hysteresis simulation in a bi-stable regime. Potential parameters K = —5, 
J = 10, L = 100, the lattice size N = 1024, and the coarsening parameter q = N. 

Table 8.1 

Approximation of the exit time r. For the statistics we use 10 4 samples and present the 95% 
confidence interval. The potential parameter J = 5, the coarse-graining parameter q = N and the 
lattice size N = 1024 are fixed. 



Parameters 


Tin 


Ttl 


Teg 


CPU m 


CPUti 


CPU cg 




microscopic 


ML-KMC 


CGMC 


[sec] 


[sec] 


[sec] 


L = N 














K = 0, /i = 1 


28.5 ±0.8 


28.3±0.8 


28.7±0.8 


1534 


9 


8 


K = 2,h = 2 


6.40±0.03 


6.40±0.03 


6.20±0.02 


884 


G 


5 


L = 100 














K = 3,h = 2.5 


6.20±0.02 


6.1±0.03 


5.93±0.02 


158 


9 


7 


K = 3,h = 3.1 


11.50±0.06 


12.4±0.1 


44.0±0.1 


526 


45 


100 



since c(x, a) — c(x, a) for all a € Sat and ieA, and c(x, a) satisfies (|4.8|) . We have 
c(x, a) — c(x, a) since 

~ f - <?{x) 

c a {x,(j) = c a {k,r))Ctf(x\k,r]) = d Q (q - rj(k)) — - 

q - 7](k) 

= d (I - <r(x)) = c a (x, a) 



and 



c d (x,<7) = c d (k, V )4(x\k,v) =rf r?(fc)e-^)4Se-^^)-^)] 

f7(fc) 
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Figure 8.4. Comparing the probability density function of the exit time from the initial coverage 
co = to c > 0.99. Potential parameters K = 2, J = 5, h = 2, L = N , the lattice size N = 1024 
and the coarse- graining parameter q = N . 



Table 8.2 

CPU time (seconds): The evolution final time T = 20, the potential parameters K = 1, J = 5, 
h = 2.5, L = N , and the coarse- graining parameter q = N 



Lattice size N Null event ML-KMC 



512 


9 


0.5 


1024 


33 


0.9 


2048 


131 


1.7 


4096 


514 


4 


8192 


2143 


13 



A. 2. Approximate dynamics. The approximate reaction rates c(x, a) denned 
in (|4.13j) satisfy the DB condition with invariant measure jiN,i3{d<j) (|4.15j) . Indeed, if 
we denote c a (x,a) = Ca(k,rj)c^ { (x\k,ri) and Cd(x,a) = c ( i{k 1 rf)c^ { (x\k, r/) we have 



c a (x,a)e-^"^^ [c a (k,r,)c^(x\r,)} e^ H ^ 



do(q - v{k)) 



1 - a(x) 



q - T}(k)_ 
d (l - <7{x))e-W {x > a) 

d a x (x)e- fir ^ x '' T) ] , 



-l3{H N (a !c )-(2a(x)-l)U{x,a)) 
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Figure 8.5. Comparing the probability density function of the exit time from the initial coverage 
co = to c > 0.99. Potential parameters K = 3, J = 5, h = 2.5, L = 100, the lattice size N = 1024, 
and the coarse- graining parameter q = N . 



and similarly 



d ri(k)e 



[do(l-o*(aj))] e -^ N(CTX) 



- / 3(// JV (o' x )-(2 -( 2; )-l)C/(x, C r)) 



Appendix B. Proof of Lemma l5.5i For the sake of completeness we also give 
the proof of Lemma 15.51 which was proved in [19 . 

Proof. We denote by V (J 0(cr) = (d x <f>(a)) xeA and c(a) = (c(x,a)) xeA . The 
equation (|5.11|) can be rewritten as 

d t u(t, a) + c(er) • V a u{t, a) = . 

For the discrete difference c^tt we obtain the equation 

d t (d x u(t, a)) + c(a) ■ V a d x u{t, a) + d x c(a) ■ V a d x u(t, a x ) = . 

From the definition of the rates c{x, a), (|4.1[) . we can estimate upper bounds for d x c{o~). 
Using that d x U(z,a) = U(z, a x ) - U(z, a) = iiC(x-^)(l-2cr(a;)) + J(x-^)(l-2cr(a;)) 



30 



E. Kalligiannaki, M. A. Katsoulakis and P. Plechac 



Two-level 

Null event 




6 8 10 12 14 16 18 20 22 24 26 

Exit time 



Figure 8.6. Comparing the probability density function of the exit time from the initial coverage 
co = to c > 0.99. Potential parameters K = 3, J = 5, h = 3.1, L = 100, the lattice size N = 1024, 
and the coarse- graining parameter q = N. Coarse-graining error of CGMC that appears due to 
the finite-range interactions is substantially reduced with the ML-KMC method (see the text for 
explanation of the inset). 



when z ^ x, and d x U(z, a) = when z = x we have 
d x c(z,(j) = < 



0(1), forz = x, 

0(1), for0< \z-x\ <S, 

0(1/ L), ioi S < \z-x\ <L. 



Then, since Cv(u) = c(a) ■ V CT w(cr), we can write 

d t (d x u(t, a)) + Cd x u{t, a) + Y d x c(z, a)d z u(t, a x ) = , 

zeA 

d t (d x u(t, a)) + C8 x u{t, a) + 0{l)d x u{t, a x ) 
+0(1) Y d z u(t,v x ) + 0(j) E d z u(t,<? x )=0- 

\z-x\<S S<\z-x\<L 

Furthermore, we have 

||9x«(t,-)lloo < \\d x u{Q, -)\\co + J 0(l)\\d x u(s,-)\\ 00 ds+ 
+ f 0(1) Y \\9xu(s,-)\\ 00 ds+ C O(j) Y l|0xu(s,-)lloods. (B.l) 

\z-x\<S 1 S<\z-x\<L 
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Figure 8.7. Average coverage trajectory. As external field h approaches critical value h c = 4 the 
coarse-graining error in CGMC becomes important. However, the two-level ML-KMC simulations 
capture correctly transitions even when CGMC fails. Potential parameters K = 3, J = 5, h = 3.1, 
L = N , the lattice size N = 1024, and the coarse- graining parameter q = N . The inset depicts the 
estimated autocorrelation function. 

Based on this relation, application of Gronwall's inequality for 6{t) = ^2 X \\d x u(t, 
and the fact that S is finite and small we conclude 

£ \\B x u(t, OIU < e c(T - t} E H^ u (°' Olloo • 

X X 

□ 

Appendix C. Kinetic Monte Carlo algorithms. Stochastic simulation algo- 
rithm (SSA) as proposed in [TT] is described next, where the evolution from a state 
er to a x is sampled by: 

Algorithm 2. Stochastic Simulation Algorithm. 

Step 1: Update, (a) Calculate all rates c(y,a),\fy £ Ajv form j4-l\ l, that are af- 
fected from the previews event. 

(b) Calculate X x (cr) = Y, y<x c (v^ a ) and A ( CT ) = E y eA« c (y^ a ) 
Step 2: Search. Obtain a uniform random number u € [0, 1) and search for x G A^r 
such that 

K-i{o-) < X(a)u < X x (a) 

Step 3 Update time from an exponential law with the parameter \(o~) or equivalently 
with the mean At = jj^- That is select a uniform random number u\ G [0, 1) 
and update the time as t' = t + St with St = — log(ui)A£. 
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The n-fold way (or BKL) algorithm, [3J, is equivalent to the SSA in the sense that 
it always leads to a successful event and requires the updating of all transition rates. 
The BKL algorithm was already designed to reduce the cost of the searching process 
by dividing the transition states into classes (the number of classes n <C N) with the 
same probability, thus the search algorithm cost depends on the number of processes 
at each site, i.e., scales linearly with the reaction range (the number of interacting 
neighbours) in the model under consideration. 

Algorithm 3. n-fold algorithm (BKL). 
Given a 

Step 1: Update, (a) Calculate all rates c(y,o~),\/y £ Ajy that are affected by the 
previous event. 

Step 2: Search. Group sites x £ A in classes D i; i = 1, . . . , n by classifying them 
with their rate values and define 

j j 
QiW = E E c(y,<7) = £|A|c(y,<r) 

i=l yeDi i=l 

for some y £ Di, j = 1, . . . , n. Generate a uniform random number u £ [0, 1) 
and search for i = 1, . . . , nsuch that 

Qi-l(cr) < Q n (cr)u < Qi(cr) , 

then choose x £ Di uniformly. 
Step 3 Update time from the exponential law with the parameter A(<r) = Q n {a), or 
equivalently with the mean At = . 

The previous two algorithms are in the class of rejection-free methods as the embeded 
Markov chain always jumps into a new state. However, by applying the uniformization 
we obtaine a null-event algorithm in which the embeded chain has nonzero probability 
to stay at the same state in each step. 

Algorithm 4. Null-event algorithm. 
Find the bounds A*' loc = do max{l, e~@ u }, and U* = min X)(r U(x, a). 
Given a 

Step 1: Search/Update. Select a site x £ A with the uniform probability i? and 
calculate c(x, a). 

Step 2:Accept/Reject. Obtain a uniform random number u £ [0,1), 

if c(x,o~) > A*' loc- u accept and update the state a — >• a x at the site x, 
ifc(x,a) < \*' loc u assign the new state to be a. 

Step 3 Update time from the exponential law with the parameter A* ,loc , or equiva- 
lently with the mean At = A „ ] loc . 
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